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' Abstract 

We investigate the surface width W of solid-on-solid surfaces in the vicin- 
ity of the roughening temperature T^. Above Tr, W"^ is expected to di- 
verge with the system size L like InL. However, close to a clean InL 
behavior can only be seen on extremely large lattices. Starting from the 
Kosterlitz-Thouless renormalization group, we derive an improved for- 
mula that describes the small L behavior on both sides of T^. For the 
Discrete Gaussian model, we used the valleys-to-mountains-reflections 
cluster algorithm in order to simulate the fluctuating solid-on-solid sur- 
face. The base plane above which the surface is defined is an L x L square 
lattice. In the simulation we took 8 < L < 256. The improved formula 
fits the numerical results very well. From the analysis, we estimate the 
roughening temperature to be Tj. = 0.755(3). 



1 Introduction 



Solid-on-solid (SOS) models are useful as interface models They belong to a large 
class of models that are believed to be in the Kosterlitz-Thouless (KT) universality 
class For SOS models, the KT transition is the roughening transition. It is 
still a challenge to devise methods for the accurate study of this transition and for 
unambiguous tests of the KT theory. 

As a prototype of an SOS model we consider the Discrete Gaussian (DGSOS) 
model, defined by the partition function 

z = T.^^p\-7^T.ih^-hA , (1) 

h [ {1,3) J 

where hi are integer "height" variables defined on the sites i of an L x L square 
lattice with periodic boundary conditions. A configuration h can be viewed as a 
surface without overhangs, embedded in three dimensions; its energy is obtained by 
summing over all nearest neighbor pairs T is the temperature (Boltzmann's 

constant is set to one). The square of the average surface width W is defined by 

W' = ^,Y.<i^h,-h,f> ■ (2) 

The model has two phases. At low temperatures the surface is smooth, and W stays 
finite as L — >■ oo. When T is increased, we encounter the roughening transition at 
T = Tj-. The KT theory predicts 0] that in the thermodynamic limit 

^ {T,-T)-^ r^\n^ (3) 

as T approaches from below is the correlation length). For T >T,., W diverges 
as L — > oo. The prediction for asymptotically large L is the "free field" behavior 
(i.e. Continuous Gaussian - hi is real instead of integer) 

^ ^ In L + const. (4) 

71 

Teff is called the "effective temperature" , and 

Teff=- for T = Tr. (5) 

TT 

Furthermore, as T approaches Tj. from above, 

Teff--~(T-Tr)5. (6) 

71 

In principle, these formulas could be used in a numerical study in order to verify 
or disprove the KT theory. In practice however this is problematic. In the smooth 
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phase, we would need unrealistically large lattices in order to test the power law 
This problem is related to the difficulties encountered in the study of the dual 
(Villain, XY) spin models, where it is hard to cleanly distinguish an essential singu- 
larity in the correlation length ^ (as predicted by KT) from a power law singularity 

^. In the rough phase, for large enough temperatures, the behavior (^ could 
be unambiguously verified numerically ||^. However, it turns out that close to Tr 
a clean logarithm is only seen on very large lattices, and in order to extract the 
values of Teg in practice we need to know the corrections to eq. (^). Otherwise we 
cannot determine Tf by checking eq. (^). Furthermore, for the largest lattice sizes 
accessible with present day computers and algorithms, it turns out that eq. (^) is 
not yet fulfilled for the region where eq. (|) holds. Actually, the status of eq. (^) is 
even worse, as will be argued later. 

In order to overcome these problems, we developed a renormalization group (RG) 
improved formula for the dependence of on L. This is our main theoretical 
result. The numerical part of our work shows that the improved formula can be 
used for extracting Tes as close as desired to Tj., from numerical simulations on 
reasonably sized lattices. We mention that very high accuracy simulations were 
possible because we have a cluster algorithm that is free of critical slowing down (the 
valleys-to- mountains- reflections algorithm 0). Vectorization |^ also helped. From 
our analysis, the best estimate for the roughening temperature is T^. = 0.755(3). 

In what follows, we first derive our improved formula, then present the analysis 
of the numerical results, and finally make some additional remarks and present our 
conclusions. 



2 RG improved finite L formula for the surface 
width 

The RG flow of the DGSOS model can be described in an x — ?/ diagram The 
trajectories are parametrized by t, the logarithm of the changing length scale. x{t) 
is related to the scale dependent ("running") temperature T{t), x(t) = TiT{t) — 2, 
while y{t) is a constant times the fugacity 0. The KT flow equations are [0: 

y{t) = -<t)yit) 

xit) = -yitf . ^ ^ 

The trajectories are hyperbolas, characterized by the constant E which depends on 
the temperature T of the model (not on the running T{t) !): 

y{tf-x{tf = E. (8) 
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Denoting e = sign{E)y\E\, and xq = x(0), the full solution of eq. (|^) is: 

E<0: x{t) = e(l+ , ^ s'^^^o I \\ 

^' V (xo + e)exp(2et) - (xo - e)y 

E = ^: x{t)= ^ (9) 

E>0: x{t) = e ^o-etan(et) 
^ ' e + Xq tan(er) 

The trajectories in the rough phase reach the free field theory, and have e < 0; in 
the smooth phase they have e > 0; at the KT transition the critical trajectory has 
e = [§. Notice that in the rough phase Tes = T(t = oo) = (2 — e)/7r. 

In order to use the RG for computing the surface width, we need to know the 
contributions corresponding to each length scale. Eq. (|^) shows that W"^ is a sum of 
a two-point-function over all distances. When increasing the lattice size L, we get 
additional additive contributions from distances of order L. Let us choose 

t = \n^, (10) 

with Lq some reference length scale, and let us approximate the sum in eq. (|^) by 
an integral. For the free field theory, the additional contributions to W"^ coming 
from an infinitesimal change in L are easily computed: since eq. (^ is always true, 
with Tcfj replaced by the temperature T, we have dW^/dt = T/n. In the case of 
the DGSOS model, the KT flow shows that for trajectories in the rough phase we 
are close to the free field theory if t is large enough. Moreover, we are also close to 
the free field theory for trajectories in the smooth phase, provided that L is much 
smaller than ^ but still large enough. Thus the main contribution to dW'^/dt will be 
similar to the free field case, the only (important) difference being that we replace 
T by the running temperature T{t): 

dt 77 

Assuming that T(t) behaves according to the KT flow for length scales larger than 
Lq, we can integrate eq. (|Tll): 



= ^ [\xit) + 2]dt + C. (12) 

The constant C contains the contributions of distances smaller than Lq. Using eq. 
(^ and eq. (H), we can express dt in terms of x and dx, after which eq. ( |T2|) reduces 
to an elementary integral. We thus obtain our final formula for W^: 



3 



which has to be used in conjunction with eqs. and i^^- 

The crucial point in the derivation of our improved formula was the replacement, 
at the appropriate stage, of the temperature T with the running temperature T{t). 
While this is a commonly used procedure in field theoretical RG arguments, it is not 
completely rigorous. A more thorough argument, based on a block-spin calculation 
within the Wilson RG framework, will be presented elsewhere 0. 

3 Simulation results 

We performed simulations of the DGSOS model for ten different values of the tem- 
perature T. At each T we considered lattice sizes of L = 8, 16, 32, 64, 128 and 256. 
Typically, we generated about 2 000 000 to 2 500 000 clusters for each temperature 
and lattice size. The expectation value of the cluster size varied in the range 0.3L^ 
to 0.35L^. The whole project required about 400 hours on one CRAY Y-MP pro- 
cessor. As will become clear from the analysis, we did not need more than half of 
our runs in order to obtain the best estimate for Tj.. However, our aim was also to 
confirm our prediction for the L-dependence of and to determine the region in 
which the improved formula is really necessary. The simulation results for W"^ are 
given in Table |l[ 

As a general rule, these data are extremely well fitted by eq. (P^, with fit 
parameters e, Xq and C. Aside from the occasional statistical fluctuation, we did 
however notice that for T < 0.755 the quality of the fits deteriorated a little. The 
important results of such fits are the value of e, which characterizes the trajectory, 
and the range of L for which the fits are good, which roughly tells us where the 
model enters the KT flow. Notice that we have to decide upon a value for Lq. The 
choice of Lq only affects the values Xq and C, as can be seen after a little algebra. 
Table ^contains the fit results for e, for all our values of T and for various fit ranges. 
Clearly, for T > 0.76 the various fit ranges give compatible results. In fact the data 
for L = 256 hardly improve things here. For T < 0.755 however, the fit results for 
e sometimes change by more than two standard deviations if we remove the data 
for L = 8. It may be that for these temperatures the KT flow is well reached only 
above L = 8. 

From Table ^ our first main numerical result strikes the eye: since e > for 
T < 0.75 and e < for T > 0.76, Tj. is between 0.75 and 0.76. This result relies 
solely on the fact that eq. (|I3|) fits the data well, and on eq. (H). 

In order to give a more precise determination of Tj., we plotted for each fit range 
e versus T, with error bars, and interpolated the two curves e + error and e — error. 
The intersection of the band thus obtained with the e = line provides an estimate 
of Tj.. In Table |^ we collected these range-dependent estimates. They are quite 
consistent with one another. Thus, taking into account the above observations 
about the quality of the fits for T < 0.755, it would not be unreasonable to quote 
as our final result the value of Tj. for the L-range 16 — 256. 
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Table 1: Simulation results for W"^. 



T 


L = 8 


L= 16 


L = 32 


L = 64 


L = 128 


L = 256 


0.740 


0.51429(34) 


0.66739(33) 


0.81589(32) 


0.96185(33) 


1.10558(34) 


1.24908(39) 


0.745 


0.51984(37) 


0.67560(34) 


0.82571(32) 


0.97408(33) 


1.12248(33) 


1.26710(40) 


0.750 


0.52537(37) 


0.68222(35) 


0.83542(35) 


0.98661(34) 


1.13615(35) 


1.28461(36) 


0.755 


0.53137(30) 


0.69052(34) 


0.84475(33) 


0.99877(33) 


1.15066(34) 


1.30165(38) 


0.760 


0.53733(37) 


0.69725(34) 


0.85444(33) 


1.00956(29) 


1.16440(28) 


1.31721(39) 


0.770 


0.54900(35) 


0.71185(32) 


0.87292(33) 


1.03164(36) 


1.18993(34) 


1.34726(37) 


0.780 


0.55902(36) 


0.72606(34) 


0.88941(33) 


1.05190(38) 


1.21432(36) 


1.37515(40) 


0.800 


0.58006(38) 


0.75240(34) 


0.92271(35) 


1.09176(40) 


1.26032(37) 


1.42846(40) 


0.820 


0.59963(37) 


0.77776(34) 


0.95355(34) 


1.12807(36) 


1.30248(38) 


1.47770(42) 


0.850 


0.62762(37) 


0.81304(37) 


0.99679(36) 


1.18076(38) 


1.36315(40) 


1.54571(43) 



Table 2: Fit results for the parameter e; the first row contains the L-range. 



T 


8-256 


16 - 256 


32 - 256 


8-128 


16 - 128 


0.740 


0.204(06) 


0.178(10) 


0.151(23) 


0.226(08) 


0.206(16) 


0.745 


0.171(07) 


0.139(14) 


0.165(21) 


0.168(11) 


0.060(57) 


0.750 


0.126(09) 


0.109(17) 


0.103(35) 


0.137(14) 


0.117(30) 


0.755 


0.030(38) 


-0.078(25) 


0.059(60) 


0.054(33) 


-0.102(33) 


0.760 


-0.124(10) 


-0.130(14) 


-0.128(27) 


-0.131(13) 


-0.151(21) 


0.770 


-0.211(06) 


-0.210(09) 


-0.221(17) 


-0.211(09) 


-0.205(17) 


0.780 


-0.277(05) 


-0.287(07) 


-0.278(14) 


-0.278(07) 


-0.300(12) 


0.800 


-0.387(04) 


-0.387(06) 


-0.388(11) 


-0.386(06) 


-0.387(10) 


0.820 


-0.478(03) 


-0.484(05) 


-0.494(09) 


-0.471(05) 


-0.474(09) 


0.850 


-0.601(03) 


-0.599(04) 


-0.590(08) 


-0.603(04) 


-0.601(07) 



Table 3: Tj. from the interpolated curves e(T). 



fit range 


8-256 


16 - 256 


32 - 256 


8-128 


16 - 128 


estimated Tr 


0.7555(25) 


0.7535(15) 


0.7550(30) 


0.7555(25) 


0.7515(55) 
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For a more conservative estimate of we plotted the values of e from the ranges 
8—256 and 16—256, together with their error bars, on the same plot. We interpolated 
the upper and lower envelopes of the error bars. From the intersection of the band 
thus obtained with the e = line we get the estimate Tr = 0.755(3). Notice that 
the best estimate in the literature 0, 71 = 0.7524(7), was obtained by a completely 
different method (matching with the critical block spin flow of the BCSOS model), 
that does not test directly any of the formulas derived from the KT theory. The best 
estimate by other authors 0, Tj. = 0.752(5) (from the analysis of the correlation 
length and susceptibility in the massive phase of the Villain model), is also consistent 
with the result presented here. 

At the beginning of section 2 we explicitly wrote down the t dependence of the 
running temperature T{t). With the numerically determined fit parameters e and 
xo, we can thus compute the flow of T(t) numerically. If we use x(t) instead of 
T(t), we can neatly plot the points {x(t),y(t)) inside the standard KT flow diagram. 
We can now do the following consistency check. The differences irAW'^/At = 
(7r/ln2) [W^^(2L) — Vr^(L)], shown in Table ^, should be discrete approximants of 
T(t), by eq. (|10|) and (pH^. Thus if we again plot the values of the points {x{t),y{t)), 
this time using the discrete approximation, we expect to obtain a similar diagram. 
We did this exercise, and indeed the two diagrams were almost identical. 

In the last column of Table ^ we show the values of T^s = (2 — e)/7r, obtained 
by again taking for each T > T^. the envelope of the error bars from the fits with 
L-ranges 8 — 256 and 16 — 256. Above Tj., if L is large enough for eq. (^ to hold, the 
running temperature stabilizes to the value Tgff. By looking at how the the results 
in the rows of Table |^ stabilize, we see that our data for W"^ enter the asymptotic 
regime of eq. (§) for 0.8 < T < 0.85 clearly, and for T = 0.78 just barely. For 
Tr < 0.77 however, this is far from being true, even at L = 256. Notice that in 
order to understand the validity region of eq. (§) we did need the values of for 
L = 256. More importantly, notice that our results show that the use of eqs. (^ 
and (^) for determining (like e.g. in [1^, |ll|) leads to a consistent underestimate. 

In order to test eq. (|^) , we fitted the values for Tes from Table ^. The fit was not 
at all good. We then allowed for a free power instead of the power |. The fit was 
now good, but the power was 0.60(4). Disregarding the point farthest away from Tj., 
T = 0.85, the situation did not improve: the power changed to 0.62(7). The fitted 
value for was 0.753(2) with the point T = 0.85 included, and 0.752(4) without 
it. While these values for Tj. are reasonable, the fact remains that the power ^ is 
not yet seen even as close to Tj. as our data in the rough phase are. Notice that this 
conclusion implies in particular that we cannot use eq. (^ in order to fit results in 
a region where the simple behavior (^ applies on lattices of still manageable size. 

As a last issue, let us remark that in the absence of a theory, one may be simply 
tempted to make some "reasonable" ansatz for the corrections to eq. (^. We tried 
to fit the data with a In In L correction (the coefficient in front of InlnL is the third 
fit parameter besides Tgfj and the constant). The fits were as good as those with eq. 
(|l3l), if not better. However, the values of Tes thus obtained were clearly wrong. It is 
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Table 4: The differences tt AW'^/At compared to T^s 



T 


8-16 


16 - 32 


32 - 64 


64 - 128 


128 - 256 


TcS 


0.740 


0.6939(21) 


0.6731(21) 


0.6615(21) 


0.6514(22) 


0.6504(24) 


T <T, 


U. ( 40 


U. IU0i)[ZO ) 


U.DoU4(^Zi ) 


n R79/I 
U.D / Z4(^Zi ) 


U.D ( ZK)[ZL ) 


U.DOODl^Zo j 


rp ^ f-p 


0.750 


0.7109(23) 


0.6944(22) 


0.6852(22) 


0.6778(22) 


0.6729(23) 


T <T, 


0.755 


0.7213(21) 


0.6990(22) 


0.6981(21) 


0.6884(21) 


0.6843(23) 




0.760 


0.7248(23) 


0.7124(21) 


0.7031(20) 


0.7018(18) 


0.6926(22) 


0.6777(48) 


0.770 


0.7381(21) 


0.7300(21) 


0.7194(22) 


0.7174(22) 


0.7131(23) 


0.7035(29) 


0.780 


0.7571(22) 


0.7404(22) 


0.7365(23) 


0.7361(24) 


0.7290(24) 


0.7267(35) 


0.800 


0.7811(23) 


0.7719(22) 


0.7662(24) 


0.7640(25) 


0.7620(25) 


0.7598(19) 


0.820 


0.8074(23) 


0.7967(22) 


0.7910(22) 


0.7905(24) 


0.7942(26) 


0.7900(22) 


0.850 


0.8404(24) 


0.8328(23) 


0.8338(24) 


0.8267(25) 


0.8274(27) 


0.8273(16) 



not difficult to understand the numerics behind this phenomenon. The main point 
is, however, to view this as another example of the danger of analyzing simulation 
results without a solid theoretical basis. 

Along the same lines, let us remark that we found a different modification of eq. 
(^ to also fit the data very well: instead of taking In L we took a power of In L (this 
power is the third fit parameter). The power that allowed for good fits very close 
to Tr never deviated from the value 1 by more than 10%. Nevertheless, the fitted 
values for T^s were again clearly wrong with this procedure. 

4 Conclusions 

We have derived a renormalization group improved formula for the finite size be- 
havior of the SOS surface width in the vicinity of the roughening transition. The 
improved formula was tested in a high accuracy simulation of the DGSOS model, 
and found to describe the data excellently. As a result of our analysis, we: 

— verified an important aspect of the KT scenario; 

— gave a precise determination of Tj.; 

— found the region in which eq. (^ cannot be used unless we consider much 
larger lattice sizes; 

— found that the region of applicability of eq. (^ is much smaller than 
previously assumed. 
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